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This  paper  describes  a  phenomenological  constitutive  model  for  ionomer  membranes  in  polymer  elec¬ 
trolyte  membrane  fuel  cells  (PEMFCs).  Unlike  the  existing  approaches  of  elasto-plastic,  viscoelastic, 
and  viscoplastic  model,  the  proposed  model  was  inspired  by  micromechanisms  of  polymer  deforma¬ 
tion.  The  constitutive  model  is  a  combination  of  the  nonlinear  visco-elastic  Bergstrom-Boyce  model 
and  hydration-temperature-dependent  empirical  equations  for  elastic  modulus  of  ionomer  membranes. 
Experiment  results  obtained  from  an  uniaxial  tension  test  for  Nation  NR-1 1 1  membrane  under  well  con¬ 
trolled  environments  were  compared  with  simulated  results  by  the  finite  element  method  (FEM)  and 
the  proposed  model  showed  fairly  good  predictive  capabilities  for  the  large  deformation  behavior  of  the 
Nation  membrane  subjected  to  the  uniaxial  loading  condition  in  a  wide  range  of  relative  humidity  and 
temperature  levels  including  liquid  water. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Mechanical  behavior  of  polymer  electrolyte  membrane  (PEM) 
has  been  studied  over  a  wide  range  of  temperature  and  relative 
humidity  (RH)  conditions  by  many  research  groups  due  to  the  sig¬ 
nificance  of  the  membrane  durability  in  fuel  cell  operation  [1-9]. 
Perfluorosulfonic  acid  (PFSA)  membrane,  originally  developed  by 
DuPont  with  the  trade  name  Nafion®,  is  widely  used  in  the  PEM 
fuel  cell.  It  is  believed  that  membrane  failure  is  the  major  fuel  cell 
life  determining  factor  [10,11].  PFSA  membrane  can  absorb  water 
and  undergoes  microstructure  transformation,  primarily  the  for¬ 
mation  of  hydrophilic  ionic  clusters  in  hydrophobic  perfluorinated 
ionomers  matrix.  This  imparts  the  membrane  proton  conductiv¬ 
ity;  it  also  causes  the  mechanical  property  of  ionomer  membrane 
to  vary  with  respect  to  temperature  and  hydration  level.  At  low 
temperature  below  90  °C,  water  acts  as  plasticizer  softening  the 
membrane  and  reducing  load  carrying  capability.  However,  at  ele¬ 
vated  temperature  above  90  °C  surprisingly,  the  opposite  trend  is 
observed;  the  more  water  a  membrane  absorbs,  the  stronger  the 
membrane  becomes  [3,5,8].  This  abnormal  behavior  was  attributed 
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to  change  in  viscoelastic  response  of  Nafion  due  to  microphase 
structural  transitions  driven  by  changes  in  temperature  and  water 
activity  [5]. 

It  is  very  important  to  know  the  maximum  level  of  membrane 
stress  when  designing  a  fuel  cell  device.  Stress  distribution  in  the 
membrane  is  also  important  for  the  understanding  of  the  mem¬ 
brane  degradation  mechanisms.  The  authors  recently  investigated 
the  mechanical  stress  effect  on  chemical  degradation  of  ionomer 
membrane  and  found  that  the  mechanical  stress  accelerates  the 
rate  chemical  decomposition  of  the  PEM  by  oxygen  radical  attacks 
[10].  Direct  measurement  of  membrane  stress  in  fuel  cells  is  not 
always  possible;  using  numerical  models  to  predict  mechanical 
stress  in  the  ionomer  membrane  in  various  conditions  in  fuel  cells 
is  often  the  only  choice.  The  macromolecular  nature  of  polymer 
is  characterized  by  the  covalently  bonded  and  long  chain  struc¬ 
ture  [  1 2  ].  The  physical  properties  of  polymeric  systems  are  strongly 
affected  by  chain  microstructure,  i.e.,  isomerism,  which  is  the  orga¬ 
nization  of  atoms  along  the  chain  as  well  as  the  chemical  identity 
of  monomer  units  [13].  An  important  feature  controlling  the  prop¬ 
erties  of  polymeric  materials  is  polymer  architecture.  The  Nafion® 
membrane  is  a  copolymer  containing  at  least  two  monomers,  i.e., 
a  TFE  backbone  and  perfluoro(4-methyl-3,  6-dioxa-7-octene-l- 
sulfonyl  fluoride)  [14].  A  large  amount  of  polymer  research  works 
continue  to  be  directed  towards  the  study  of  molecular  mechanisms 
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Fig.  1.  Tensile  stress-strain  curve  of  Nafion  NRE  212  tested  at  room  temperature 
with  4.2  mm  s-1  pulling  rate. 


governing  their  structure-property  relationships.  Among  them,  the 
stress-strain  response  of  polymers  has  been  recognized  for  a  long 
time  as  one  of  the  most  informative  properties  [15].  Fig.  1  shows 
the  typical  stress-strain  curve  for  a  Nafion  NRE-212  membrane. 
Macroscopic  nature  of  the  mechanical  behavior  for  the  Nafion® 
membrane  under  the  tensile  stress  before  rupture  is  character¬ 
ized  by  an  elastic  response  (Hook’s  law),  followed  by  the  strain 
hardening  in  the  plastic  deformation  range  after  the  yield  point. 
These  elastic  and  plastic  deformation  for  the  membrane  is  also 
time-dependent,  i.e.,  visco-elastic  and  visco-plastic.  Experimental 
data  presented  in  the  Solasi’s  work  [16]  clearly  demonstrated  the 
complicated  non-linear  time-,  hydration-level-,  and  temperature- 
dependent  behavior  of  the  ionomer  membrane. 

It  is  assumed  that  when  an  external  load  is  applied  to  a  polymer, 
the  molecular  bonds  experience  stress,  and  in  order  to  relieve  them¬ 
selves  as  much  as  possible,  the  chain  segments  undergo  internal 
rearrangements  [15];  the  way  the  polymer  reacts  to  the  exter¬ 
nal  stress  is  dependent  on  the  magnitude  and  rate  of  the  applied 
stress,  chain  morphology,  environmental  factors  such  as  relative 
humidity  (RH)  and  temperature.  In  literatures,  it  is  believed  that  the 
Nafion®  membrane  consists  of  at  least  two  phases  [17];  an  amor¬ 
phous  phase  and  a  crystalline  phase,  and  the  crystallinity  for  1100 
EW  membrane  is  in  a  range  between  5  and  20%  [18].  Therefore,  it 
is  expected  that  each  component  contributes  to  the  deformation 
resistance  differently. 

Haward  and  Thackray  [19]  attempted  to  interpret  this 
macroscopic  constitutive  behavior  of  polymer  based  on  the 
microstructure  of  polymers.  The  polymer’s  mechanical  response 
was  described  by  with  two  parallel  processes  as  the  initial 
non-linear  elastic  response  controlled  by  the  secondary  inter- 
molecular  interactions,  and  the  entanglement  network  response 
controlled  by  the  primary  intramolecular  interactions  which 
give  an  entropic  contribution  at  large  strains.  As  a  continuous 
attempt  for  describing  the  mechanical  behavior  of  polymers, 
specifically,  ionomer  membrane,  we  propose  a  nonlinear  visco¬ 
elastic-visco-plastic  constitutive  model  by  introducing  hydration 
dependence  into  Bergstrom  and  Boyce’s  model  [20]  published  in 
1998. 

2.  Constitutive  modeling  of  ionomer  membrane 

2.1.  Micro-mechanism  of  polymer  deformation 

The  constitutive  model  of  ionomer  membranes  is  needed  for 
continuum  mechanics  model  to  predict  the  distribution  of  the 


Amorphous  region 


Crystallite 


Fig.  2.  Schematic  of  hierarchy  of  ionomer  membrane  structure 


stress  and  strain  in  the  fuel  cell  membrane.  The  typical  stress-strain 
behavior  in  Fig.  1  can  be  qualitatively  described  by  the  initial 
Hookean  elasticity  until  the  stress  developed  becomes  sufficiently 
large  to  produce  a  plastic  deformation  at  the  imposed  rate.  Follow¬ 
ing  a  strain  hardening  region  after  the  yield  point,  the  deformation 
eventually  leads  to  the  final  rupture  of  the  membrane;  it  is  gen¬ 
erally  believed  that  the  rupture  of  the  polymer  is  induced  by  a 
defect  (microcracks,  crazing)  formation  and  accumulation  process 
[21,22].  To  model  the  mechanical  behavior  of  the  polymer  mem¬ 
brane  subjected  to  uniaxial  stretching,  researchers  have  applied 
linear  elasticity  [23],  visco-plasticity  [24],  and  elasto-plastic  the¬ 
ory  [1,6,25].  None  of  the  models  can  account  for  the  effects  of  all 
known  micro-mechanisms  of  the  ionomer  deformation,  such  as  a 
reptational  plastic  flow,  chain  entanglement,  and  entropic  effect  on 
deformation. 

Fig.  2  shows  a  hypothetical  hierarchical  structure  of  the 
ionomer  membrane  consisting  of  the  amorphous  network,  crys¬ 
talline  regions,  and  water  clusters.  It  has  been  established  by  many 
investigators  that  molecular  chain  orientation  evolves  with  mag¬ 
nitude  and  state  of  strain  to  produce  strain  hardening  in  polymers 
[26,27]  and  the  polymer  macromolecular  structure,  which  forms  a 
network  by  virtue  of  physically  entangled  molecular  chains,  evolves 
during  deformation  as  secondary  valence  interactions  dissociate 
with  plastic  strain  [28].  Entanglements,  a  topological  constraint, 
develop  from  the  interpenetration  of  random-coil  chains  and  are 
of  importance  in  determining  rheological,  dynamic,  and  fracture 
properties  [29,30]. 

Chain  mobility  and  morphological  relaxation  of  intermolecu- 
lar  chain  are  expected  to  increase  with  temperature,  which  can  be 
explained  by  reptational  dynamics  [31  ].  Also,  it  is  reported  that  the 
mechanical  stress  increases  the  molecular  mobility  during  plastic 
deformation  [32]. 

Considering  the  characteristic  morphology  and  microstructure 
of  the  ionomer  membrane  reported  by  a  number  of  researches 
[15,33,34],  the  mechanical  deformation  mechanisms  of  the  PFSA 
based  ionomer  are  hypothesized  as  follows: 

1.  The  physical  crosslinks  formed  by  molecular  entanglements, 
ionic  interaction  in  sulfonic  acid  groups  and  intermolecular,  sec¬ 
ondary  (Van  der  Waals)  interactions  in  crystalline  phase  bears 
the  low  stresses  in  the  elastic  region  of  the  stress-strain  curve, 
leaving  the  major  portion  of  covalent  bonds  unaffected. 

2.  When  the  external  load  is  increased  beyond  a  certain  level, 
ionic  domains  start  to  permanently  deform,  elongate  and  reori¬ 
ent  [34];  as  the  stress  develop,  some  chains  in  the  amorphous 
and/or  crystalline  phase  can  overcome  the  secondary  inter¬ 
actions  and  develop  irreversible  slippage  and  reorientation, 
yielding  occurs.  As  stress/strain  increases  further,  the  perma- 
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Fig.  3.  One  dimensional  rheological  representation  of  the  constitutive  model  for  the 
ionomer  membrane. 


nently  entangled  chains  (which  can  not  slip  out  of  physical 
entanglements)  become  taut  and  start  locking  up,  resulting  in 
strain  hardening.  As  strain  increases  further,  small  crystallites 
can  disintegrate  [33]. 

Based  on  the  hypothesized  deformation  micro-mechanisms 
for  the  ionomer  membrane,  a  ID  rheological  representation  of 
the  constitutive  model  (as  shown  in  Fig.  3)  is  proposed.  The 
mechanical  behavior  for  ionomer  membranes  can  be  decom¬ 
posed  into  two  parts:  a  visco-plastic  response,  which  is  relevant 
to  irreversible  molecular  chain  slippage  and  a  time-dependent 
viscoelastic  response.  The  viscoelastic  response  can  be  further 
decomposed  into  the  response  of  two  molecular  networks  acting  in 
parallel:  a  first  network  (A)  represents  the  equilibrium  of  the  visco¬ 
elastic  response  and  a  second  network  (B)  the  time-dependent 
deviation  from  the  visco-elastic  equilibrium  state.  This  idea,  called 
Dual  Network  Fluoropolymer  (DNF)  model,  was  initially  developed 
by  Bergstrom  for  modeling  fluoropolymers  [35].  The  decompo¬ 
sition  idea  was  also  used  by  Boyce  [36]  and  Bergstrom  [37,38]. 
The  Cauchy  stress  acting  on  the  network  A  and  B  can  be  modeled 
by  any  of  the  classical  nonlinear  hyper-elasticity  models  for  elas¬ 
tomers.  In  this  research,  the  Cauchy  stress  which  is  a  function  of  the 
Cauchy-Green  deformation  tensor  followed  that  of  the  Bergstrom 
and  Boyce’s  model  [20]  for  elastomers,  which  is  based  on  the  eight- 
chain  model  of  Arruda  and  Boyce  [26].  Also,  the  plastic  flow  rule  for 
the  network  B  is  motivated  by  reptational  dynamics  of  a  polymer 
[31,39].  The  DNF  model  does  not  account  for  the  hydration  effect 
on  the  mechanical  properties  of  ionomer  membranes  originated 
from  the  water  containing  ionic  clusters  formed  by  sulfonic  acid 
groups.  For  the  modeling  of  ionomer  membranes,  it  is  assumed  that 
hydration  effect  can  be  incorporated  implicitly  into  the  empirical 
equation  for  the  elastic  modulus,  which  is  explained  in  details  in 
the  modeling  section  below. 


2.2.  Constitutive  modeling 

For  the  analysis  of  the  large  deformation  of  polymers,  the  con¬ 
cept  of  multiplicative  decomposition  of  the  deformation  gradient 
into  elastic  and  plastic  parts  had  been  typically  employed  instead 
of  an  additive  decomposition  [27,36-38,40-44].  The  mathemati¬ 
cal  description  of  the  constitutive  model  for  ionomer  membranes 
is  based  on  the  breakdown  of  the  overall  deformation  into  the 
visco-elastic  and  visco-plastic  deformation  which  is  referred  to  as 
the  Kroner-Lee  decomposition.  The  deformation  gradient  Fean  be 
expressed  by  Eq.  (1),  as  shown  below  [35]: 

F  =  FveFp  (1) 

F?  is  the  deformation  due  to  pure  plastic  flow  representing 
irreversible  chain  motion,  and  Fve  =  FA  =  FB  is  the  remaining  con¬ 
tribution  to  F  associated  with  distortion  and  reorientation  of 
crystallites  and  entanglement  of  polymer  chains. 

The  visco-elastic  deformation  gradient  is  further  decomposed 
into  elastic  and  viscous  parts: 

Fve=FeFv  (2) 

Here,  F  is  the  reversible  (elastic)  deformation  gradient  and  Fv 
indicates  the  viscous  deformation  gradient.  The  spatial  velocity  gra¬ 
dient  L  is  given  by  L  =  F  •  F_1 .  By  inserting  F  =  FveFp  into  L  =  F  •  F-1 , 
the  corresponding  rate  kinematics  can  be  decomposed  into  visco¬ 
elastic  and  visco-plastic  contributions  as  well: 

L  =  F  ■  F-1  =  ( FveFP  +  FveFP)(FveFP)~ 1  =  Lve  +  V>  (3) 

Here,  Lp  =  Dp  +  Wp  can  be  described  as  the  sum  of  the  symmetric 
(the  rate  of  deformation,  Dp )  and  the  skew-symmetric  part  (the  spin 
tensors,  Wp.)  Similarly,  the  velocity  gradient  of  viscoelastic  parts 
can  be  decomposed  into  elastic  and  viscous  components:  Lve  = 
Fve .  Fve~x  =  Le  +  Ly,  where  Lv  =  Dv  +  Wv.  The  intermediate  config¬ 
urations  described  by  p  and  v ,  in  general,  cannot  be  uniquely 
determined,  since  an  arbitrary  rigid  rotation  can  be  superimposed 
on  it  and  leave  it  stress  free  [45].  The  intermediate  state  can  be 
determined  uniquely  in  different  ways  and  one  convenient  way  is 
to  prescribe  Wv  =  0  and  Wp  =  0,  which  means  that  the  flow  is  irro- 
tational  [46].  In  addition  to  that,  plastic  and  viscous  deformation 
are  assumed  to  be  incompressible,  i.e.,  det(Fy)  =  1  and  det(FP)  =  l. 
In  our  study,  the  volumetric  swelling  and  shrinkage  behavior  of  the 
Nation  as  functions  of  the  hydration  level  and  temperature  are  not 
considered  in  the  kinematics  during  the  deformation,  the  hydration 
level  and  temperature  are  assumed  to  be  constants  in  a  given  load¬ 
ing  condition.  Also,  the  deformation  of  Nation  over  whole  strain 
range  is  assumed  to  be  nearly  incompressible,  det(F)^  1  as  well. 

The  Cauchy  stress  tensor  for  network  A  is  given  by  the  eight- 
chain  representation  [26,37]: 

TA=f8ch(Fve)=-^=  ■  ViX?dev  [B"e* *]  +*[/*-  1 11  W 

jve\ve  L-1(l/A/oc/c) 

where  ]ve  =  det  [Fye],  fiA  is  a  temperature  and  hydration  level 
dependent  initial  shear  modulus,  Xlock  is  the  chain  locking  stretch, 
Bve*  =  (jvey2/3p/e(pvey  j s  th< e  left  Cauchy  Green  tensor,  Xve  = 

Bn/3  is  the  effective  chain  stretch  based  on  the  eight- 
chain  assumption,  L_1(x)  is  the  inverse  Langevin  function,  where 
L(x)  =  coth  (x)  -  1  /x,  and  k  is  the  bulk  modulus.  To  obtain  the  inverse, 
a  curve  fit  of  the  inverse  Langevin  function  is  used  for  all  x  as  follows 
[39]: 

U  f  1. 31446  tan(l. 58986* )  +  0.91209*,  if  |x|  <  0.84136  ^ 

*  f  l/(sign(x)-x)  ifO.84136  <  |x|  <  1 

The  Cauchy  stress  tensor  for  network  B  can  be  calculated  from  the 
same  eight-chain  representation  that  was  used  for  network  A  and 
computed  by  multiplication  the  eight-chain  expression  on  the  elas¬ 
tic  deformation  gradient  F  with  a  scalar  factor  Sb  which  can  be 
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Table  1 

Fitting  constants  for  elastic  modulus  equation. 


Al 

Bx 

a2 

b2 

0.000645 

-0.058673 

-0.014673 

10.534189 

considered  as  a  specific  material  parameter. 


TB=Sb  ■ f8ch(Fe)=SB  ■ 


L~\XeIXlock) 


dev[Be*]  +  K\]e 


Then,  the  total  Cauchy  stress  can  be  calculated  as  the  sum  of  the 
two  resultants,  T=TA  +  TB.  The  First  Piolar  Kirchhoff  stress,  P  which 
relates  forces  in  the  present  configuration  with  areas  in  the  ref¬ 
erence  configuration  can  be  calculated  from  conversion  between 
two  stress,  P=J ■  T ■  F~T  =J  •  {TA  +  TB)  •  F~T .  In  order  to  incorporate  the 
hydration-  and  temperature-dependent  mechanical  property  such 
as  nA  and  for  ionomer  membranes,  empirical  relationship  of 
elastic  modulus  as  a  function  of  membrane  water  content,  Am  and 
temperature,  0  (°C)  was  used.  The  bulk  modulus  k  is  assumed  to  be 
a  constant  and  its  magnitude  was  taken  as  a  reasonable  value  in  the 
analysis.  It  was  found  that  the  simulated  stress  is  not  significantly 
affected  by  the  bulk  modulus.  Earlier  attempt  for  the  empirical  rela¬ 
tion  of  the  elastic  modulus  was  made  by  Hsu  et  al.  [47].  In  [48],  the 
uniaxial  testing  data  for  N1 1 1  membrane  were  collected  under  con¬ 
trolled  temperature  and  RH  environments  and  these  data  sets  were 
used  for  fitting  the  exponential  type  function  as  shown  below: 


F(A,  0)  =  exp  {(At  •  0  +  )  •  Am  +  (A2  ■  0  +  B2 )} 


(7) 


As  reported  in  [50],  another  way  to  eliminate  the  numerical  dif¬ 
ficulty  is  to  use  certain  differentiable  smooth  ramp  function  R(a) 

— V  — V 

(shown  in  Fig.  4.)  which  replaces  [A  -  1  ]  by  R(A  -  1 ) 
where, 

R(a)= 2.8853-e-  [  0.3466  +  0.1733-  -  0.5  •  In  i  ,  ^  „  *  ,  ,  „  1  )  (12) 

v  ;  y  e  \  cosh[0.3446  ■(«/£)] //  v  J 

satisfying 


R{  0)  =  £, 


=  1. 


where  Am  =  0.043  +  ]7.8]aT  -  39.85a^  +  36. 0a^  forO  <  aT  <  1 
[49]  and  aT  is  the  water  activity  (RH)  defined  by  aT  =  PwlPsat{0 ), 
where  pw  is  water  vapor  pressure  and  psat  is  saturation  water  vapor 
pressure  at  the  temperature.  From  the  elastic  modulus  equation, 
the  initial  shear  modulus  for  network  A  and  B  can  be  described  by 


As  a  result,  the  velocity  gradient  of  the  viscous  flow  can  be 
expressed  as: 

Lv  =  FeLvFe~'  =Dv  =  yv  — 
re 


HB  =Sob  -E(X,0)  =  Sob  ■  exp{(+  ■  6  +  B^)  ■  Xm +(A2  ■  6  +  B2)}  (8)  _i  (  t*  \ 

F  =F  I  Y  ]  F  F 

Ma  =  s0a  -  Mb  (9)  V  x  / 


(13) 


where  s0a  and  s0B  are  material  parameters  and  Alf  f?i,  A2,  and  B2 
are  fitting  constants  listed  in  Table  1. 

The  rate  of  visco-plastic  flow  of  network  B  can  be  described 
by  Dv  =  yvNv.  The  tensor  Nv  specifies  the  directions  of  the  driv¬ 
ing  stresses  of  the  relaxed  configuration  convected  to  the  current 
configuration,  and  the  terms  yv  indicates  the  flow  rates  being  given 
by  the  reptation-inspired  equation  [37]: 


Yv  =  Yo^V-  1 


( _ — _ 

V  T- base  +  fiPe 


n 


(10) 


where  re  =  ||Te  ||  =  (tr[Te  Te  ])  is  the  Frobenius  norm  of  Te  = 
dev[Te],  the  direction  of  the  driving  stress  is  described  by  Nv  = 
Te  /re,  Xv  =  y //tr(By*)/3  is  an  effective  viscous  chain  stretch, 
gv*  _  yvy2/3pv(pvjT  js  ^h e  left  Cauchy-Green  deformation  ten¬ 
sor,  pe  =  -(T^  +  T|2  +  T|3)/3  is  the  hydrostatic  pressure,  and 
Cg  [-1,0],  m  >  0,  n,  yo,  6base  and  rbase  are  material 

— V  C 

parameters.  As  pointed  out  by  Bergstrom,  the  term  [A  -  1  ]  cap¬ 
tures  a  strain-dependence  of  the  effective  viscosity  and  this  might 
cause  the  term  to  grow  numerically  very  large  if  the  effective  stretch 
A  is  unity  or  close  to  unity  in  both  the  unloaded  state  and  when 
the  applied  strain  switches  between  tension  and  compression.  One 
way  to  resolve  this  problem  is  to  introduce  a  parameter  £  &  0.01  to 
avoid  the  singularity  and  Eq.  (10)  can  be  modified  as  follows: 


y"  =  Yo[2’  -  1  +  e]C  ■  ( 


T base  +  fiPe 


n 


(11) 


The  rate  of  plastic  flow  can  be  described  by  a  phenomenological 
equation  [35]: 

•  p  =  f  ab(e-e0)b_1e  if  z  >  cr0 

Y  |  0  otherwise, 

where,  a  >  0,  b  >  0  and  cr0  >  0  are  material  parameters,  r  =  1 1  dev[T]  1 1  f 
is  the  Frobenius  norm  of  the  deviatoric  part  of  the  Cauchy  stress  T, 
and  £0  is  the  effective  strain  when  r  is  equal  to  <j0;  the  effective 
strain  can  be  calculated  from  £  =  1 1  Eln  II,  where  Fln  =  ln[V]  and  V 
is  the  left  stretch  tensor,  and  £  is  the  effective  strain  rate.  For  our 
research,  £0  is  considered  as  a  constant  and  the  engineering  strain 
rate  was  used  for  £  for  the  simplicity  since  the  plastic  flow  rate  can 
be  controlled  by  choosing  the  appropriate  parameters,  a  and  b.  As 
can  be  noticed  in  Eq.  (14),  the  plastic  flow  rate  is  a  function  of  the 
strain  rate  and  the  magnitude  of  current  strain. 

In  summary,  the  velocity  gradient  of  the  plastic  flow  can  be 
expressed  as: 

IP  =  pve^ppve-'1  =  QP  =  yP  dev[T] 

r 

FP  =  Fve~'  (yPde^T1)  FveFp  (15) 

3.  Experimental 

The  experimental  data  of  Nation  that  we  used  to  fit  and  verify  the 
proposed  model  came  from  data  reported  in  Yue  Zou’s  thesis  work 
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Fig.  5.  Schematic  diagram  of  membrane  mechanical  testing  setup. 


[48].  The  mechanical  test  was  performed  on  MTS  Tytron™  250,  a 
horizontal  load  frame  designed  especially  for  membrane  testing. 
An  environment  chamber  (T  and  RH)  was  built  to  fit  the  horizon¬ 
tal  rail  of  Tytron  load  frame  and  the  RH  control  was  achieved  by 
mixing  dry  and  saturated  gas  stream.  Steady-state  RH  values  were 
taken  from  a  chilled  mirror  dew  point  sensor  (EdgeTech  Dew  Prime 
II)  continuously  sampling  the  RH  of  the  gas  existing  the  chamber 
and  temperature  was  measured  by  a  platinum  RTD  probe  installed 
closely  to  the  membrane  sample  inside  the  chamber.  All  the  uni¬ 
axial  tension  tests  data  used  in  this  study  were  conducted  under 
environment  with  well-controlled  relative  humidity  and  temper¬ 
ature.  The  membrane  samples,  Nation  NR-111  (hereby  referred 
to  as  Nlll  sample),  were  cut  into  rectangular  shape  with  about 
50  mm  length  and  6.5  mm  wide.  The  specimens  were  equilibrated 
in  the  chamber  for  at  least  1  h  before  the  tests.  Mechanical  test¬ 
ing  data  at  three  different  condition,  which  were  25  °C  and  80%RH, 
25  °C  and  50%RH,  and  65  °C  and  75%RH  were  used  for  compari¬ 
son  with  FEM  results.  The  stress-strain  curve  for  the  samples  is 
obtained  by  applying  a  tensile  force  at  a  uniform  strain  rate  of 
0.0132  s-1. 

To  study  the  material  properties  in  liquid  water  hydrated  state,  a 
water  bath  tray  together  with  a  U-type  pulling  rod  were  designed  to 
perform  tensile  testing  with  the  membrane  sample  fully  immersed 
in  water.  The  nominal  cross-sectional  area  and  gauge  length  of 
the  specimens  were  used  in  the  calculation  of  the  engineering 
stress-strain  curves  from  the  displacement-load  data  for  the  vapor- 
equilibrated  membrane  and  a  linear  expansion  rate  of  15%  [51] 
were  applied  to  compensate  the  volumetric  expansion  in  the 
calculation  of  the  cross-section  area  for  the  water-equilibrated 
membrane  at  80  °C.  The  engineering  stress  and  strain  are  converted 
to  true  stress  and  true  strain  when  compared  with  the  model  pre¬ 
diction  (Fig.  5). 

4.  Finite  element  simulation 

The  proposed  constitutive  model  for  ionomer  membrane  was 
implemented  into  Comsol  Mutiphysics  3.5  software  package.  The 
Structure  Mechanics  and  the  PDE  modules  were  utilized  for  the 
simulation  of  time  dependent  behavior  of  the  membrane  in  an 


application  module  of  plane  stress.  The  PDE  module  was  used  for 
the  integration  of  the  time  evolution  equation  for  viscous  flow  and 
plastic  flow.  At  every  time  step,  the  viscous  and  plastic  deforma¬ 
tion  gradients  were  calculated  and  used  for  computing  the  Cauchy 
stress  acting  on  the  network  A  and  B.  To  model  the  nearly  incom¬ 
pressible  material  for  ionomer  membranes,  mixed  U-P  formulation 
provided  in  the  software  was  used  for  calculation  of  an  independent 
variable,  pressure  p. 

5.  Results  and  discussion 

Fig.  6(a)  shows  the  engineering  stress-strain  curves  of  as 
received  Nlll  membrane  tested  at  three  different  environmental 
conditions  and  strain  rates.  Plastic  deformation  sets  in  at  a  strain 
around  0.1  and  the  material  hardens  along  the  rest  of  stress-strain 
curve.  At  large  strain,  the  elastic  component  of  the  strain  is  negli¬ 
gible  compared  to  the  plastic  strain.  It  can  be  seen  that  the  elastic 
modulus  and  yield  stress  are  decreasing  with  increasing  temper¬ 
ature  and  RH.  The  slope  of  the  strain-hardening  portion  of  the 
stress-strain  curve  changes  slightly  with  respect  to  temperature 
and  humidity. 

The  stress-strain  curves  from  FEM  simulation  (Fig.  6b-d) 
shows  fairly  good  qualitative  agreement  with  experimental  behav¬ 
ior  of  vapor-equilibrated  ionomer  membranes.  The  FEM  predicts 
temperature-  and  hydration-dependent  mechanical  behavior  of 
membrane  and  the  overall  shape  of  the  strain-hardening  behavior 
matches  the  experiment  results.  The  proposed  constitutive  model 
can  accurately  describe  the  mechanical  behavior  of  the  vapor- 
equilibrated  membrane.  Material  parameters  and  other  constants 
used  at  the  simulation  are  tabulated  in  Tables  2  and  3. 

From  the  rheological  representation  in  Fig.  3,  the  true  stress 
components  can  be  dissociated  into  the  stress  acting  on  network  A 
and  network  B.  As  mentioned  earlier,  the  true  stress  from  network 
A  captures  the  equilibrium  response  of  the  material  and  network 
B  represents  time-dependent  deviation  from  viscoelastic  equilib¬ 
rium  state.  As  generally  accepted,  physical  resistances  which  are 
related  with  intermolecular  and  intramolecular  interactions  gov¬ 
ern  the  energy  barrier  that  must  be  overcome  to  yield  the  material 
and  to  deform  it  up  to  large  plastic  strain  [52].  These  molecular 
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Engineering  strain 


Engineering  strain 


Fig.  6.  (a)  Experiment  results  of  stress-strain  curves  of  Nafion  N1 1 1  membrane 
25  °C  and  80%RH,  and  (d)  65  °C  and  75%RH. 


Engineering  strain 


various  test  conditions,  (b)  FEM  simulation  results  of  N1 1 1  at  25  °C  and  50%RH,  (c) 


Table  2 

Material  parameters  used  in  the  FEM  simulation  of  vapor-equilibrated  membranes. 


c 

-0.5 

m 

6 

n 

4.5 

P 

0.6 

Yo 

1 

a 

0.01 

b 

0.78 

so 

0.01 

cr0 

5  MPa 

SB 

8 

e 

0.0132 

@base 

100 

K 

97  MPa 

hlock 

6 

interactions  induce  the  rate  and  temperature  effects  prevailing  in 
the  material  behavior. 

Figs.  7  and  8  show  the  true  stress  components  from  the 
stress-strain  curve  calculated  from  FEM.  At  low  strain,  the  stress 
is  dominantly  exerted  by  the  intermolecular  interactions  (network 
B)  in  crystallites,  amorphous  phase,  and  ionic  domain.  However, 


Table  3 

T  and  RH  dependent  material  parameters  used  for  FEM  simulation  of  vapor- 
equilibrated  membranes. 


25  °C,  50%RH 

25  °C,  80%RH 

65  °C,  75%RH 

Those 

Mb*  1.42 

Mb*  1.42 

Mb*  3 

So A 

0.75 

0.8 

1 

Sob 

0.0256 

0.0384 

0.033 

at  high  strains,  network  B  can  no  longer  bear  the  stress  due  to 
the  molecular  relaxation  and/or  crystallite  disintegration,  and  the 
rubber-like  network  A  starts  to  dominate.  This  results  in  the  strain 
hardening  behavior,  which  can  be  explained  by  deformation,  reori¬ 
entation,  and  tightening  of  entangled  chain  molecules.  In  Fig.  7,  the 
results  showed  that  at  the  same  temperature  but  different  RH  (25  °C 
50%  and  80%),  the  mechanical  behavior  of  the  rubber-like  network 
A  does  not  show  discernable  difference,  but  the  whole  curve  for 
network  B  is  shifted  downward  as  the  RH  increases,  which  indi¬ 
cates  that  the  water  vapor  weakens  the  network  B  component  more 


Fig.  7.  Comparison  of  contribution  of  each  stress  components  for  viscoelastic  net¬ 
work  A  and  B  from  FEM  simulation  results  for  N1 1 1  at  25  °C,  50%  and  80%RH. 
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Fig.  8.  Contribution  of  each  stress  components  of  viscoelastic  network  A  and  B  from 
FEM  simulation  results  for  N1 1 1  at  25  °C,  80%RH  and  65  °C,  75%RH. 


than  network  A.  This,  in  turn,  imply  that  the  water  vapor  interferes 
with  intermolecular  network  such  as  secondary  interaction  of  PTFE 
backbone  and  strong  ionic  interaction  in  the  hydrophilic  domain, 
and  deteriorates  their  interactions  [53  ].  However,  as  shown  in  Fig.  8, 
the  temperature  appears  to  influences  critically  on  both  mechan¬ 
ical  behavior  of  network  A  and  B  at  the  same  RH  and  reduces 
load-carrying  capability  of  membranes  by  softening  intermolecu¬ 
lar  interactions  and  enhancing  the  chain  mobility  in  the  material. 
There  has  been  recent  attempt  to  interpret  the  mechanical  behavior 
by  microstructure  transition  [5,9]  due  to  the  temperature  and  water 
activity,  and  it  is  reported  that  the  combined  effects  of  temperature 
and  water  alter  the  structure  of  the  hydrophilic  domains  chang¬ 
ing  the  number,  strength  and  flexibility  of  cross-links  between 
domains.  The  attractive  interactions  between  sulfonic  acid  groups 
are  likely  to  induce  ionic  group  aggregation,  which  form  effective 
cross-links  that  stiffen  membrane  at  the  low  temperature.  How¬ 
ever,  increase  of  temperature  causes  the  sulfonic  acid  groups  to 
become  randomly  dispersed  and  thus  breaks  these  cross-links.  The 
mechanical  behavior  predicted  by  our  FEM  analysis  is  consistent 
with  these  hypotheses. 

In  order  to  evaluate  the  effectiveness  of  the  proposed  constitu¬ 
tive  model,  water-equilibrated  Nation  samples  at  80  °C  were  tested 
and  the  stress-strain  curves  from  the  uniaxial  tests  were  used  to  fit 
the  model.  In  addition  to  that,  two  different  strain  rates  of  0.3  s-1 
and  0.0045  s-1  were  used  to  verify  that  the  model  can  capture 
the  rate-dependent  mechanical  behavior  of  the  ionomer  mem¬ 
brane.  Elastic  modulus  of  water-equilibrated  membrane,  E,  was 
determined  from  the  linear  curve  fit  to  the  experimental  results 
measured  from  the  uniaxial  tension  test.  The  material  parameters 
are  listed  in  Tables  4  and  5,  and  the  results  of  the  stress-strain 


Table  4 

Material  parameters  used  in  the  FEM  simulation  of  water-equilibrated  membranes 
at  80  °C. 


c 

-0.5 

m 

6 

n 

4.5 

P 

0.6 

Yo 

1 

b 

0.78 

So 

0.01 

S0A 

2.1 

SB 

4 

@base 

100 

K 

97  MPa 

^lock 

6 

E 

25  MPa 

Table  5 

Strain-rate  dependent  material  parameters  used  in  the  FEM  simulation  of  water- 
equilibrated  membranes  at  80  °C. 


Strain  rate,  0.3 

Strain  rate,  0.0045 

'Zbase 

liB*  2 

IjLb*  2.52 

a 

0.1 

0.2 

&o 

5  MPa 

3  MPa 

Sob 

0.077 

0.0625 

Fig.  9.  Comparison  the  experimental  data  from  uniaxial  tension  test  with  FEM 
results  (a)  at  the  strain  rate  0.3  s-1  and  (b)  the  strain  rate  0.0045  s_1  under  the  water 
at  80  °C. 


curves  are  plotted  in  Fig.  9.  The  simulated  stress-strain  curves  show 
that  the  model  can  predict  the  rate-dependence  of  the  mechanical 
behavior  as  expected.  The  true  stress  components  of  network  A 
and  B  at  the  two  different  strain  rate  are  plotted  in  Fig.  10,  showing 
that  high  strain  rate  stiffens  both  the  time-dependent  network  A 
and  B. 

Finally,  the  stress  acting  on  each  component  A  and  B  are  mainly 
determined  by  the  magnitude  of  initial  shear  modulus  fiA  which 
is  hypothetically  related  to  the  rubber-like  network  such  as  the 
molecular  chain  entanglement  in  amorphous  phase  of  ionomer 
membrane  and  is  responsible  for  the  hardening  behavior,  and  /x# 
being  associated  with  intermolecular  interaction  including  ionic 
clusters,  polymer  backbone,  and  side  chains.  It  is  observed  that 
Sqa,  a  ratio  of  / iA  and  /xB  is  less  than  unity  for  the  FEM  analysis 


Fig.  10.  Stress  components  from  viscoelastic  network  A  and  B  as  revealed  by  the 
FEM  simulation  results  for  N1 1 1  at  80  °C  under  the  water  at  strain  rate  of  0.3  s-1  and 
0.0045  s-1. 
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of  the  vapor-equilibrated  membrane  indicating  that  there  is  initial 
stiff  response  before  the  viscous  flow  effects  dominates  the  behav¬ 
ior.  The  magnitude  of  fiA  is  found  to  be  almost  two  times  greater 
than  that  of  /X#  for  the  analysis  of  membrane  immersed  in  a  liquid 
water.  Furthermore,  the  true  strain  at  which  the  true  stress  compo¬ 
nents  from  network  A  and  B  intersect  is  approximately  0.7  for  the 
vapor  equilibrated  membranes  as  can  be  seen  in  Figs.  7  and  8.  How¬ 
ever,  the  intersection  point  for  the  water  equilibrated  membranes 
is  shifted  down  to  a  value  less  than  0.3.  These  observations  sug¬ 
gest  that  liquid  water  considerably  deteriorates  interactions  in  the 
hydrophilic  domains  and  any  secondary  interactions  which  might 
be  present  in  the  crystalline  and  amorphous  phase  being  respon¬ 
sible  for  initial  stiffness,  and  mechanical  behavior  of  the  water 
equilibrated  membrane  exhibits  more  likely  the  behavior  of  elas¬ 
tomeric  materials,  hence  rubber-like  behavior  start  dominates  the 
stress-strain  response  at  a  much  smaller  strain  than  that  of  the 
vapor-equilibrated  membranes. 


6.  Conclusion 

We  have  developed  a  new  constitutive  model  describing  the 
finite  deformation  of  the  ionomer  membrane  for  PEMFCs.  The 
constitutive  relationship  is  nonlinear  viscoelastic-viscoplastic, 
and  strain-rate,  temperature,  and  hydration  dependent. 
Micromechanism-inspired  Bergstrom-Boyce  model  was  employed 
to  capture  the  nonlinear  viscoelastic  behavior  of  the  membrane 
and  reptational  dynamics  inspired  flow  rule  for  viscous  flow  was 
used  to  describe  the  time  dependent  plastic  behavior.  The  pro¬ 
posed  model  can  describe  the  stress-strain  behavior  of  both  vapor- 
and  liquid-water-equilibrated  membranes,  and  rate-dependent 
mechanical  behavior.  The  simulated  stress  from  FEM  analysis  is 
obtained  from  summation  of  two  different  molecular  networks 
acting  in  parallel.  Network  A  produces  strain-hardening/stiffening 
behavior  resulting  from  molecular  reorientation,  entanglement, 
and  locking  up  due  to  the  large  deformation  and  the  network 
B  generates  the  initially  stiff  response  as  well  as  the  rate,  tem¬ 
perature,  and  hydration  dependence  of  initial  flow.  It  was  found 
that  water  softens  the  network  B  component  and  tempera¬ 
ture  affects  significantly  the  material  behavior  of  both  network 
components. 

By  choosing  the  appropriate  material  parameters,  this  model 
can  accurately  capture  the  mechanical  behavior  of  ionomer  over 
a  wide  range  of  temperature  and  hydration  level,  implying  that 
the  hypothesized  mechanisms  govern  the  material  behavior  of 
ionomer  membrane.  Future  work  is  needed  to  improve  the  consti¬ 
tutive  model  by  incorporating  the  volumetric  expansion  as  function 
of  hydration  level  into  the  kinematic  equation,  and  predicting 
the  stress  and  strain  responses  of  membranes  subjected  to  com¬ 
plex  loading  conditions  and  history,  such  as  multi-axial  stress  in 
constrained  membrane  during  dry  out.  As  indicated  earlier,  we 
have  assumed  the  water  content  and  temperature  in  the  mem¬ 
brane  remain  constant  in  the  deformation  process.  The  authors 
realize  that  this  is  an  idealized  situation.  In  reality,  due  to  the 
non-linear  interactions  between  the  membrane  and  water,  the 
membrane  may  not  come  to  complete  equilibrium  with  water 
during  the  deformation  process.  In  a  fuel  cell,  the  ionomer  mem¬ 
brane  is  constantly  adjusting  itself  (deformation  is  one  aspect  of 
such  adjustment)  to  equilibrate  with  typically  non-uniform  and 
time-varying  liquid-vapor  mixture  of  water.  It  takes  finite  time 
for  water  to  be  transported  in  or  out  of  the  membrane,  in  tran¬ 
sient  operation  conditions  membrane  may  not  be  able  to  come  into 
fully  equilibrated  state  during  the  deformation  process.  Coupled 
water  transport  and  mechanical  deformation  models  are  necessary 
to  accurately  the  time-dependent  water-stress-strain  state  in  the 
ionomers  membrane. 
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